# ============ Section 0: load packages & data ============ 
library(dplyr)
library(tidyr)
library(psych)
library(ggplot2)
library(cregg)
library(miceadds)
library(sensemakr)
library(jtools)
library(sjPlot)
library(mediation)

load("IncivilityDoesNotExist_data.RData")

# ============ Section 1: Pre-processing: data cleaning, reshaping, computing variables ============

## ============ Section 1.1: Data cleaning ============

# Remove speeders, straightliners, and those who failed attention checks

DAT.CLEAN.INTERIM <- IncivilityDoesNotExist_data %>% 
  dplyr::mutate(
    # Compute the median value of `DURATION`
    DURATION_MEDIAN = median(DURATION), 
    # Create a new variable `SPEEDERS`: 'Yes' if respondent completed survey in less than half of the median duration, otherwise 'No'
    SPEEDERS = if_else(DURATION < as.numeric(DURATION_MEDIAN)/2, "Yes", "No")
    ) %>% 
  dplyr::rowwise() %>% 
  dplyr::mutate(
    # Compute the standard deviation of responses to battery of questions with reversed items (authoritarian aggression)
    POLATT_AUTHORITY_SD = sd(c(POLATT_AUTHORITY_1, POLATT_AUTHORITY_2, 
                               POLATT_AUTHORITY_3, POLATT_AUTHORITY_4, 
                               POLATT_AUTHORITY_5, POLATT_AUTHORITY_6)),
    # Create  variable `STRAIGHTLINERS`: 'Yes' if the respondent gave identical answers to all authority questions (SD = 0), otherwise 'no'
    STRAIGHTLINERS = if_else(POLATT_AUTHORITY_SD == 0, "Yes", "No")
    ) %>% 
  dplyr::ungroup() 

DAT.CLEAN.INTERIM %>% 
  group_by(ATTENTIONCHECK) %>% 
  dplyr::summarise(n=n()) #number of failed attention cheks
DAT.CLEAN.INTERIM %>% 
  group_by(STRAIGHTLINERS) %>% 
  dplyr::summarise(n=n()) #number of straigthliners
DAT.CLEAN.INTERIM %>% 
  group_by(SPEEDERS) %>% 
  dplyr::summarise(n=n()) #number of speeders

# Filter out inattentive respondents respondents
DAT.CLEAN <- DAT.CLEAN.INTERIM %>% #filter out
  dplyr::filter(ATTENTIONCHECK == "Passed") %>% #failed the attention check
  dplyr::filter(STRAIGHTLINERS == "No") %>% #straight-liners
  dplyr::filter(SPEEDERS == "No") #speeders

## ============ Section 1.2: Reshaping the dataset for analyses (wide to long format)  ============

# This section reshapes the dataset from wide to long format to facilitate analysis.
# Each respondent has evaluated 6 profiles (3 tasks with 2 profiles per task).
# The dataset is restructured so that each respondent appears 6 times, once for each profile they evaluated.

# Create vectors for the variable names by combining attribute prefixes and suffixes.
# The prefixes denote different attributes and dependent variable, while suffixes indicate tasks and profiles.
prefixes <- c("ATT_PARTY", "ATT_GENDR", "ATT_TRAIT", "ATT_POPUL", "ATT_ATTAK", "ATT_ISSUE", "ATT_ISSUEALL", "ATT_SETTG", "ATT_IMPOL")
suffixes <- c("T1_P1", "T1_P2", "T2_P1", "T2_P2", "T3_P1", "T3_P2")
# Function to create full variable names by combining prefixes with suffixes
create_vector <- function(prefix, suffixes) {
  paste0(prefix, "_", suffixes)
}
# Generate vectors of variable names for each prefix
vector_party <- create_vector("ATT_PARTY", suffixes)
vector_gendr <- create_vector("ATT_GENDR", suffixes)
vector_trait <- create_vector("ATT_TRAIT", suffixes)
vector_popul <- create_vector("ATT_POPUL", suffixes)
vector_attak <- create_vector("ATT_ATTAK", suffixes)
vector_issue <- create_vector("ATT_ISSUE", suffixes)
vector_issueAll <- create_vector("ATT_ISSUEALL", suffixes)
vector_settg <- create_vector("ATT_SETTG", suffixes)
vector_impol <- create_vector("ATT_IMPOL", suffixes)
vector_polite <- create_vector("DV_IMPOLITE", suffixes)
vector_accept <- create_vector("DV_UNACCEPT", suffixes)
vector_normal <- create_vector("DV_NOTNORMAL", suffixes)
vector_right <- create_vector("DV_WRONG", suffixes)
vector_feeling <- create_vector("DV_FEELING", suffixes)

# Pivot the dataset from wide to long format for each attribute and dependent variable.
# This process transforms the dataset so that each row represents a single respondent's evaluation of a profile.
DAT.CLEAN.PVT <- cbind(
  # Pivot the data for each attribute and dependent variable
  DAT.CLEAN %>% pivot_longer(all_of(vector_party), names_to = "TASK", values_to = "ATT_PARTY"),
  DAT.CLEAN %>% pivot_longer(all_of(vector_gendr), names_to = "TASK", values_to = "ATT_GENDR") %>% dplyr::select(ATT_GENDR),
  DAT.CLEAN %>% pivot_longer(all_of(vector_trait), names_to = "TASK", values_to = "ATT_TRAIT") %>% dplyr::select(ATT_TRAIT),
  DAT.CLEAN %>% pivot_longer(all_of(vector_popul), names_to = "TASK", values_to = "ATT_POPUL") %>% dplyr::select(ATT_POPUL),
  DAT.CLEAN %>% pivot_longer(all_of(vector_attak), names_to = "TASK", values_to = "ATT_ATTAK") %>% dplyr::select(ATT_ATTAK),
  DAT.CLEAN %>% pivot_longer(all_of(vector_issue), names_to = "TASK", values_to = "ATT_ISSUE") %>% dplyr::select(ATT_ISSUE),
  DAT.CLEAN %>% pivot_longer(all_of(vector_issueAll), names_to = "TASK", values_to = "ATT_ISSUALL") %>% dplyr::select(ATT_ISSUALL),
  DAT.CLEAN %>% pivot_longer(all_of(vector_settg), names_to = "TASK", values_to = "ATT_SETTG") %>% dplyr::select(ATT_SETTG),
  DAT.CLEAN %>% pivot_longer(all_of(vector_impol), names_to = "TASK", values_to = "ATT_IMPOL") %>% dplyr::select(ATT_IMPOL),
  DAT.CLEAN %>% pivot_longer(all_of(vector_polite), names_to = "TASK", values_to = "DV_IMPOLITE") %>% dplyr::select(DV_IMPOLITE),
  DAT.CLEAN %>% pivot_longer(all_of(vector_accept), names_to = "TASK", values_to = "DV_UNACCEPT") %>% dplyr::select(DV_UNACCEPT),
  DAT.CLEAN %>% pivot_longer(all_of(vector_normal), names_to = "TASK", values_to = "DV_NOTNORMAL") %>% dplyr::select(DV_NOTNORMAL),
  DAT.CLEAN %>% pivot_longer(all_of(vector_right), names_to = "TASK", values_to = "DV_WRONG") %>% dplyr::select(DV_WRONG),
  DAT.CLEAN %>% pivot_longer(all_of(vector_feeling), names_to = "TASK", values_to = "DV_FEELING") %>% dplyr::select(DV_FEELING)) %>% 
  # Update the `TASK` variable to more descriptive names using one of the 
  mutate(TASK = case_when(TASK == "ATT_PARTY_T1_P1" ~  "T1A", 
                          TASK == "ATT_PARTY_T1_P2" ~  "T1B", 
                          TASK == "ATT_PARTY_T2_P1" ~  "T2A", 
                          TASK == "ATT_PARTY_T2_P2" ~  "T2B", 
                          TASK == "ATT_PARTY_T3_P1" ~  "T3A", 
                          TASK == "ATT_PARTY_T3_P2" ~  "T3B")) %>% 
  # Drop the original wide-format columns since they are no longer needed
  dplyr::select(-c(ATT_GENDR_T1_P1:ATT_SETTG_T3_P2, DV_IMPOLITE_T1_P1:DV_FEELING_T3_P2))  

## ============ Section 1.2: Compute variables ============

# Creation of variables based on attribute and respondents characteristics
# Computation of scales for the main dependent variable (i.e., incivility perceptions) and measured independent variables (e.g., populist attitudes)

DAT.ANALYSES <- DAT.CLEAN.PVT %>%
  
  # Compute the `ISSUE_IMPORTANCE` variable by mapping the discussed issue to the respondent's self-reported importance.
  # The `ISSUE_IMPORTANCE` levels are recoded into descriptive labels and converted to an ordered factor.
  dplyr::mutate(
    ISSUE_IMPORTANCE = case_when(
      ATT_ISSUALL == "Abortion" ~ ISSUEIMP_ABORTION, 
      ATT_ISSUALL == "IndividualRights&Liberties" ~ ISSUEIMP_RIGTSLIB, 
      ATT_ISSUALL == "GunOwnership" ~ ISSUEIMP_GUNOWNER,
      ATT_ISSUALL == "theEconomy" ~ ISSUEIMP_ECONNOMY, 
      ATT_ISSUALL == "TaxReform" ~ ISSUEIMP_TAXATION, 
      ATT_ISSUALL == "Inflation" ~ ISSUEIMP_INFLATIN),
    ISSUE_IMPORTANCE = case_when(
      ISSUE_IMPORTANCE == 1 ~ "Not at all important",
      ISSUE_IMPORTANCE == 2 ~ "Slightly important",
      ISSUE_IMPORTANCE == 3 ~ "Important",
      ISSUE_IMPORTANCE == 4 ~ "Fairly important",
      ISSUE_IMPORTANCE == 5 ~ "Very important"),
    ISSUE_IMPORTANCE = factor(ISSUE_IMPORTANCE, levels = c("Not at all important", "Slightly important", "Important", "Fairly important", "Very important")),
    
    # Create the `ATT_PARTYDIR` variable indicating whether the party in the profile matches the respondent's party.
    # "Inparty" if they match, "Outparty" if they differ, and "Unknown" if the party is unknown or missing.
    ATT_PARTYDIR = case_when(
      PID_DIR == as.character(ATT_PARTY) ~ "Inparty", 
      PID_DIR != as.character(ATT_PARTY) & as.character(ATT_PARTY) != "Unknown" ~ "Outparty",
      as.character(ATT_PARTY) == "Unknown" | is.na(PID_DIR) ~  "Unknown"),
    
    # Setting the order of other attribute levels
    ATT_GENDR = factor(ATT_GENDR, levels = c("Female", "Male")),
    ATT_TRAIT = factor(ATT_TRAIT, levels = c("High communion, low agency", "High agency, low communion")),
    ATT_POPUL = factor(ATT_POPUL, levels = c("Insider", "Outsider")),
    ATT_PARTY = factor(ATT_PARTY, levels = c("Democrat", "Republican", "Unknown")),
    ATT_PARTYDIR = factor(ATT_PARTYDIR, levels = c("Unknown", "Outparty", "Inparty")),
    ATT_ATTAK = factor(ATT_ATTAK, levels = c("Attack", "Counterattack")),
    ATT_ISSUE = factor(ATT_ISSUE, levels = c("Economic issue", "Moral issue")),
    ATT_ISSUALL = factor(ATT_ISSUALL, levels = c("theEconomy", "Inflation", "TaxReform", "GunOwnership", "IndividualRights&Liberties", "Abortion")),
    ATT_IMPOL = factor(ATT_IMPOL, levels = c("Insult", "Ridicule")),
    ATT_SETTG = factor(ATT_SETTG, levels = c("Low incivility context", "High incivility context"))
    ) %>% 
  
  # Reverse-code items for certain scales to ensure consistency in measurement.
  dplyr::mutate(
    POLATT_POPULISM_3_REV = 8 - POLATT_POPULISM_3,
    POLATT_AUTHORITY_3_REV = 8 - POLATT_AUTHORITY_3,
    POLATT_AUTHORITY_4_REV = 8 - POLATT_AUTHORITY_4, 
    POLATT_AUTHORITY_5_REV = 8 - POLATT_AUTHORITY_5) %>% 
  rowwise() %>% 
  
  # Compute mean values for dependent and independent variables by aggregating related items into scales.
  dplyr::mutate(
    
    # Compute the scale for incivility perceptions (`DV_NORMATIVE`) as the mean of related items.
    # Then compute the overall incivility scale (`DV_UNCIVIL`) by including the normative perceptions and impoliteness measures.
    DV_NORMATIVE = mean(c(DV_UNACCEPT, DV_NOTNORMAL, DV_WRONG)), 
    DV_UNCIVIL = mean(c(DV_NORMATIVE, DV_IMPOLITE)),
    
    # Compute scales for various independent variables by averaging related items.
    POLATT_POPULISM = mean(c(
      POLATT_POPULISM_1, POLATT_POPULISM_2, POLATT_POPULISM_3_REV, 
      POLATT_POPULISM_4, POLATT_POPULISM_5, POLATT_POPULISM_6, 
      POLATT_POPULISM_7
    )),
    
    POLATT_AUTHORITY = mean(c(
      POLATT_AUTHORITY_1, POLATT_AUTHORITY_2, POLATT_AUTHORITY_3_REV, 
      POLATT_AUTHORITY_4_REV, POLATT_AUTHORITY_5_REV, POLATT_AUTHORITY_6
    )),
    
    SOCEXP_MARGINALZ = mean(c(
      SOCEXP_STATSGAIN, SOCEXP_STATSLOSS
    )),
    
    PERSONL_NARRCISM = mean(c(
      PERSONL_NARRCISM_1, PERSONL_NARRCISM_2, PERSONL_NARRCISM_3, 
      PERSONL_NARRCISM_4
    )),
    
    PERSONL_PSYCHOPT = mean(c(
      PERSONL_PSYCHOPT_1, PERSONL_PSYCHOPT_2, PERSONL_PSYCHOPT_3, 
      PERSONL_PSYCHOPT_4
    )),
    
    PERSONL_MACCHIVL = mean(c(
      PERSONL_MACCHIVL_1, PERSONL_MACCHIVL_2, PERSONL_MACCHIVL_3, 
      PERSONL_MACCHIVL_4
    )),
    
    PERSONL_AGGRESSV = mean(c(
      PERSONL_AGGRESSV_1, PERSONL_AGGRESSV_2, PERSONL_AGGRESSV_3, 
      PERSONL_AGGRESSV_4, PERSONL_AGGRESSV_5, PERSONL_AGGRESSV_6, 
      PERSONL_AGGRESSV_7, PERSONL_AGGRESSV_8, PERSONL_AGGRESSV_9, 
      PERSONL_AGGRESSV_10, PERSONL_AGGRESSV_11, PERSONL_AGGRESSV_12
    ))
  ) %>% 
  # Ungroup data to return it to a non-rowwise state.
  dplyr::ungroup() 

# ============ Section 2: Sample & variables descriptive ============

# Sample description is reported in Appendix F
# Summary statitistics for dependent and independent variables are in manuscript (methods)

Hmisc::describe(
  DAT.ANALYSES %>% 
    dplyr::select(
      SEX, 
      AGE, 
      EDU_CAT1, 
      POLINT, 
      LR_SELF, 
      LC_SELF, 
      PID_DIR,
      DV_IMPOLITE, #impolite perceptions
      DV_UNACCEPT, #unacceptabile perceptions
      DV_NOTNORMAL, #not normal perceptions
      DV_WRONG, #wrong perceptions
      DV_UNCIVIL, #incivilitz perceptions
      DV_FEELING, #candidate evaluations
      POLATT_POPULISM, #populist attitudes
      POLATT_AUTHORITY, #authoritarian aggression
      SOCEXP_STATSLOSS, #perceptions of status loss
      SOCEXP_STATSGAIN, #perceptions of status gain
      SOCEXP_MARGINALZ, #perceptions of marginalization
      PERSONL_PSYCHOPT, #psychopathy
      PERSONL_MACCHIVL, #Machiavellianism
      PERSONL_NARRCISM, #narcissism
      PERSONL_AGGRESSV #trait aggression
      )
  )
    
# ============ Section 3: Hypotheses Testing (Model A1) ============

# This section conducts main hypotheses testing using the `fun_amce` function for Model A1.
# The results are presented in the manuscript (pp. 19-20, Figure 2) and in the appendix (Appendix E, Table E1).

# Define the function `fun_amce` to perform Average Marginal Component Effects (AMCEs) and marginal means estimation.
# The function uses the `cj` function from the `cregg` package to handle conjoint analysis.
fun_amce <- function(dat, fit, test){
  cj(dat,
     fit, # The formula for the model, e.g., DV_UNCIVIL ~ ATT_GENDR + ATT_TRAIT + ...
     alpha = 0.008, # Setting the significance level to 0.008 to adjust for multiple comparisons.
     id = ~ RESPID, # Clustering by respondent (RESPID) to account for repeated measures.
     weights = NULL, # No weights are applied here.
     estimate = test, # Specifies whether to estimate AMCE or marginal means ('amce' or 'mm').
     feature_order = c("ATT_GENDR", "ATT_TRAIT", "ATT_POPUL", "ATT_ATTAK", "ATT_ISSUE", "ATT_SETTG", "ATT_PARTYDIR"),  # Specifies the order of the features in the output.
     feature_labels = list( # Assigning hypothesis labels for each attribute (e.g., "H1a" for ATT_GENDR).
       ATT_GENDR = "H1a", 
       ATT_TRAIT = "H1b", 
       ATT_POPUL = "H1c", 
       ATT_ATTAK = "H2a", 
       ATT_ISSUE = "H2b", 
       ATT_SETTG = "H3", 
       ATT_PARTYDIR = "Party ID"
     )
  )
}

# Fit models to estimate both AMCEs and marginal means using `fun_amce`.
# These models are used to generate results for the main hypotheses testing.

# Estimate Average Marginal Component Effects (AMCEs) for the dependent variable `DV_UNCIVIL`
MA1_AM_FIT <- fun_amce(DAT.ANALYSES, 
                       DV_UNCIVIL ~ ATT_GENDR + ATT_TRAIT + ATT_POPUL + 
                         ATT_ATTAK + ATT_ISSUE + ATT_SETTG + ATT_PARTYDIR, 
                       "amce")
# Estimate marginal means for the same model
MA1_MM_FIT <- fun_amce(DAT.ANALYSES, 
                       DV_UNCIVIL ~ ATT_GENDR + ATT_TRAIT + ATT_POPUL + 
                         ATT_ATTAK + ATT_ISSUE + ATT_SETTG + ATT_PARTYDIR, 
                       "mm")

# Plot the results of AMCE estimates
# This plot is included in the manuscript (Figure 2)
MA1_AMCE_PLOT <- 
  plot(MA1_AM_FIT, feature_headers = FALSE) + 
  facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "left") + 
  theme(legend.position = "none") +
  scale_color_manual(values = c("black", "black", "black", "black", "black", "black", "black"))

# Plot the results of marginal means
# This plot is also be included in the manuscript (Figure 2)
MA1_MM_PLOT <- 
  plot(MA1_MM_FIT, feature_headers = FALSE, vline = 55) + 
  facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "left") + 
  xlim (50,63) +
  theme(legend.position = "none") +
  scale_color_manual(values = c("black", "black", "black", "black", "black", "black", "black"))

# Combine AMCE and marginal means plots into a single figure for comparison
# This combined plot is included in the manuscript as Figure 2
egg::ggarrange(MA1_AMCE_PLOT + theme(legend.position = "none"), 
               MA1_MM_PLOT + theme(legend.position = "none", axis.text.y = element_blank()), 
               nrow = 1)

## ============ Section 3.1: Robustness check 1; OLS regression model (Model A2, Appendix E, Table E2) ============

# This section performs the first robustness check by replicating the main analysis (MA1) using OLS regression models.
# The results of this robustness check are presented in Appendix E, Table E2.

# Fit OLS regression model (Model A2) using clustered standard errors.
# The model estimates the effect of various attributes on the dependent variable `DV_UNCIVIL`.
# Clustering by respondent ID (`RESPID`) accounts for repeated measures.
MA2_OLS_FIT <- miceadds::lm.cluster(
  DV_UNCIVIL ~ 
    ATT_GENDR + ATT_TRAIT + ATT_POPUL + ATT_ATTAK + ATT_ISSUE + ATT_SETTG + 
    ATT_PARTYDIR + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE),
  cluster = "RESPID", 
  data = DAT.ANALYSES)

# Export the the OLS regression results to a formatted table for the appendix.
# This table is included in Appendix E, Table E2.
export_summs(MA2_OLS_FIT$lm_res,
             coefs= c("Male (Candidate)" = "ATT_GENDRMale", 
                      "High agency,low communion" = "ATT_TRAITHigh agency, low communion",
                      "Outsider" = "ATT_POPULOutsider", 
                      "Counterattack" = "ATT_ATTAKCounterattack", 
                      "Moral issue" = "ATT_ISSUEMoral issue", 
                      "High incivility context" = "ATT_SETTGHigh incivility context", 
                      "Outparty" = "ATT_PARTYDIROutparty",
                      "Inparty" = "ATT_PARTYDIRInparty",
                      "Male (Respondent)" = "SEXMale", 
                      "Political interest" = "as.numeric(POLINT)", 
                      "Party ID" = "as.numeric(PID_STR)",
                      "Issue importance" = "as.numeric(ISSUE_IMPORTANCE)",
                      "(Intercept)"),
             ci_level = 0.99,
             stars = c(`***` = 0.001, `**` = 0.005, `*` = 0.008))

## ============ Section 3.2: Robustness check 2; Incivility type (Model A3 & A4, Appendix E, Figure E1) ============

# This section performs the second robustness check by replicating the main analysis (MA1) for different types of incivility.
# Results are presented in Appendix E, Figure E1.

# Fit models to estimate Average Marginal Component Effects (AMCEs) by incivility type (ATT_IMPOL: insult vs. ridicule).
# This model calculates the AMCEs for various attributes, splitting results by the type of incivility.
MA34_AM_FIT <- cj(
  DAT.ANALYSES, DV_UNCIVIL ~ ATT_GENDR + ATT_TRAIT + ATT_POPUL + ATT_ATTAK + 
    ATT_ISSUE + ATT_SETTG + ATT_PARTYDIR, 
  id = ~RESPID, 
  estimate = "amce", 
  by = ~ ATT_IMPOL, # Splitting the results by the type of incivility (insult vs. ridicule).
  feature_order = c("ATT_GENDR", "ATT_TRAIT", "ATT_POPUL", "ATT_ATTAK", "ATT_ISSUE", "ATT_SETTG", "ATT_PARTYDIR"),
  feature_labels = list(ATT_GENDR = "H1a", ATT_TRAIT = "H1b", ATT_POPUL = "H1c", ATT_ATTAK = "H2a", ATT_ISSUE = "H2b", ATT_SETTG = "H3", ATT_PARTYDIR = "Party ID"))

# Fit models to estimate Marginal Means (MMs) by incivility type (ATT_IMPOL: insult vs. ridicule).
# This model calculates the marginal means for various attributes, splitting results by the type of incivility.
MA34_MM_FIT <- cj(
  DAT.ANALYSES, DV_UNCIVIL ~ ATT_GENDR + ATT_TRAIT + ATT_POPUL + ATT_ATTAK + 
    ATT_ISSUE + ATT_SETTG + ATT_PARTYDIR, 
  id = ~RESPID, 
  estimate = "mm", 
  by = ~ATT_IMPOL,
  feature_order = c("ATT_GENDR", "ATT_TRAIT", "ATT_POPUL", "ATT_ATTAK", 
                    "ATT_ISSUE", "ATT_SETTG", "ATT_PARTYDIR"),
  feature_labels = list(ATT_GENDR = "H1a", ATT_TRAIT = "H1b", ATT_POPUL = "H1c", 
                        ATT_ATTAK = "H2a", ATT_ISSUE = "H2b", ATT_SETTG = "H3", 
                        ATT_PARTYDIR = "Party ID"))

# Plot the AMCEs by incivility type (Appendix E, Figure E1).
# This plot displays how the AMCEs differ by the type of incivility (insult vs. ridicule) for each attribute.
MA34_AM_PLOT <- plot(
  MA34_AM_FIT, feature_headers = FALSE, legend_title = "Incivility type") + 
  facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "left") + 
  theme(legend.position = "right") +
  aes(shape = BY) +
  scale_shape_manual(name = "Incivility type", values=c(1, 2)) + 
  scale_colour_manual(values=rep("black", 7))

# Plot the MMs by incivility type (Appendix E, Figure E1).
# This plot displays how the marginal means differ by the type of incivility (insult vs. ridicule) for each attribute.
MA34_MM_PLOT <- plot(
  MA34_MM_FIT, feature_headers = FALSE, legend_title = "Incivility type", vline = 50) + 
  facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "left") + 
  theme(legend.position = "right")  +
  aes(shape = BY) +
  scale_shape_manual(name = "Incivility type", values=c(1, 2)) + 
  scale_colour_manual(values=rep("black", 7))

# Combine the AMCE and MM plots into a single figure for easier comparison (Appendix E, Figure E1).
# The combined figure shows both AMCEs and MMs for each attribute split by incivility type.
egg::ggarrange(MA34_AM_PLOT + 
                 theme(legend.position = "bottom", legend.justification = "left"), 
               MA34_MM_PLOT + 
                 theme(legend.position = "none", axis.text.y = element_blank()), 
               nrow = 1) 

## ============ Section 3.3: Robustness check 3; Impoliteness & inappropriateness as DVs (Models A5, A6, A7, A8, Appendix E, Figure E2) ============

# This section performs a robustness check by replicating the main analysis (MA1) for various indicators of incivility perceptions.
# The indicators include:
# - Model A5: "impolite"
# - Model A6: "wrong"
# - Model A7: "unacceptable"
# - Model A8: "not normal"
# Results are presented in Appendix E, Figure E2.

# Fit AMCEs models for each indicator of incivility perceptions using the custom function `fun_amce`.
# This function estimates the Average Marginal Component Effects (AMCEs) for different dependent variables.
MA5_AMCE_FIT <- fun_amce(DAT.ANALYSES,
                         DV_IMPOLITE ~ ATT_GENDR + ATT_TRAIT + ATT_POPUL + 
                           ATT_ATTAK + ATT_ISSUE + ATT_SETTG + ATT_PARTYDIR, 
                         "amce")
MA6_AMCE_FIT <- fun_amce(DAT.ANALYSES, 
                         DV_UNACCEPT ~ ATT_GENDR + ATT_TRAIT + ATT_POPUL + 
                           ATT_ATTAK + ATT_ISSUE + ATT_SETTG + ATT_PARTYDIR, 
                         "amce")
MA7_AMCE_FIT <- fun_amce(DAT.ANALYSES, 
                         DV_WRONG ~ ATT_GENDR + ATT_TRAIT + ATT_POPUL + 
                           ATT_ATTAK + ATT_ISSUE + ATT_SETTG + ATT_PARTYDIR, 
                         "amce")
MA8_AMCE_FIT <- fun_amce(DAT.ANALYSES, 
                         DV_NOTNORMAL ~ ATT_GENDR + ATT_TRAIT + ATT_POPUL + 
                           ATT_ATTAK + ATT_ISSUE + ATT_SETTG + ATT_PARTYDIR, 
                         "amce")

# Fit MMs models for each indicator of incivility perceptions using the "mm" option in `fun_amce`.
# This function estimates the Marginal Means (MMs) for different dependent variables.
MA5_MM_FIT <- fun_amce(DAT.ANALYSES, 
                       DV_IMPOLITE ~ ATT_GENDR + ATT_TRAIT + ATT_POPUL + 
                         ATT_ATTAK + ATT_ISSUE + ATT_SETTG + ATT_PARTYDIR, 
                       "mm")
MA6_MM_FIT <- fun_amce(DAT.ANALYSES, 
                       DV_UNACCEPT ~ ATT_GENDR + ATT_TRAIT + ATT_POPUL + 
                         ATT_ATTAK + ATT_ISSUE + ATT_SETTG + ATT_PARTYDIR, 
                       "mm")
MA7_MM_FIT <- fun_amce(DAT.ANALYSES, 
                       DV_WRONG ~ ATT_GENDR + ATT_TRAIT + ATT_POPUL + 
                         ATT_ATTAK + ATT_ISSUE + ATT_SETTG + ATT_PARTYDIR, 
                       "mm")
MA8_MM_FIT <- fun_amce(DAT.ANALYSES, 
                       DV_NOTNORMAL ~ ATT_GENDR + ATT_TRAIT + ATT_POPUL + 
                         ATT_ATTAK + ATT_ISSUE + ATT_SETTG + ATT_PARTYDIR, 
                       "mm")

# Combine AMCEs results into a single dataframe for further processing and visualization.
MDVS_AMCE_FIT <- rbind(MA5_AMCE_FIT, MA6_AMCE_FIT, MA7_AMCE_FIT, MA8_AMCE_FIT) %>% 
  # Recoding outcome variable names for readability.
  mutate(outcome = case_when(
    outcome == "DV_IMPOLITE"  ~ "Impolite",
    outcome == "DV_UNACCEPT" ~ "Unacceptable",
    outcome == "DV_WRONG" ~ "Wrong",
    outcome == "DV_NOTNORMAL" ~ "Not normal"),
    # Setting the order of the outcome variable for consistent plotting.
    outcome = factor(outcome, levels = c("Impolite", "Wrong", "Unacceptable", "Not normal"))) %>% 
  # Filtering to remove reference levels which are always set to zero.
  filter(level %in% c("Male","High agency, low communion","Outsider", 
                      "Counterattack", "Moral Issue", "High incivility context",
                      "Inparty", "Outparty"))

# Combine MMs results into a single dataframe for further processing and visualization.
MDVS_MM_FIT <- rbind(MA5_MM_FIT, MA6_MM_FIT, MA7_MM_FIT, MA8_MM_FIT) %>% 
  # Recoding the outcome variable names, similar to AMCE model
  mutate(outcome = case_when(
    outcome == "DV_IMPOLITE"  ~ "Impolite",
    outcome == "DV_UNACCEPT" ~ "Unacceptable",
    outcome == "DV_WRONG" ~ "Wrong",
    outcome == "DV_NOTNORMAL" ~ "Not normal"),
    # Setting the order of the outcome variable for consistent plotting.
    outcome = factor(outcome, levels = c("Impolite", "Wrong", "Unacceptable", "Not normal"))) 

# Plotting AMCEs results for each indicator of incivility perceptions.
# This plot displays how the AMCEs vary by the type of dependent variable (impolite, wrong, unacceptable, not normal).
MDVS_AM_PLOT <- plot(
  MDVS_AMCE_FIT, feature_headers = FALSE, group = "outcome", legend_title = "Dependent variable") + 
  facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "left") +
  aes(shape = outcome) +
  scale_shape_manual(name = "Dependent variable", values=c(1, 2, 3, 4)) + 
  scale_colour_manual(values=rep("black", 7))

# Plotting MMs results for each indicator of incivility perceptions.
# This plot displays how the marginal means vary by the type of dependent variable (impolite, wrong, unacceptable, not normal).
MDVS_MM_PLOT <- plot(
  MDVS_MM_FIT, feature_headers = FALSE, group = "outcome", vline = 50, legend_title = "Dependent variable") + 
  facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "left") +
  aes(shape = outcome) +
  scale_shape_manual(name = "Dependent variable", values=c(1, 2, 3, 4)) + 
  scale_colour_manual(values=rep("black", 7))

# Combine the AMCE and MM plots into a single figure for easy comparison (Appendix E, Figure E2).
# The combined figure allows for side-by-side comparison of AMCEs and MMs for each indicator of incivility perceptions.
egg::ggarrange(MDVS_AM_PLOT + 
                 theme(legend.position = "bottom", legend.justification = "left"), 
               MDVS_MM_PLOT + 
                 theme(legend.position = "none"), 
               nrow = 1) 

# ============ Section 4: Exploratory Analysis; Individual disposition (Model B1) ============

# This section conducts the exploratory analyses on individual differences for Model B1.
# The results are presented in the manuscript (pp. 20-22, Figure 3 & 4) and in the appendix (Appendix E, Table E3).

# Define function to create forest plots for model results
fun_forest <- function(fit, labels, colors, models_labels, title, title2){
  plot_models(fit, 
              rm.terms = c("AGE", "SEXMale", "EDU_NUM", "as.numeric(POLINT)", "as.numeric(PID_STR)", "as.numeric(ISSUE_IMPORTANCE)", "ATT_PARTYDIROutparty", "ATT_PARTYDIRInparty"),
              axis.labels = labels,
              axis.title = title,
              axis.lim = c(-3,2),
              title = title2,
              legend.title = "Model",
              m.labels = models_labels,
              show.values = TRUE, 
              value.size = 3,
              dot.size = 2,
              show.p = TRUE,
              vline.color = "grey",
              colors = colors) +
    theme_bw() +
    theme(
      plot.title = element_text(size = 10, face = "bold", hjust = 0.50),
      axis.text.y = element_text(hjust = (0)),
      axis.text.x = element_text(size = 10),
      text = element_text(size = 10, face = "italic"), 
      legend.position = "right",
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank())
}

# Define labels and colors for the forest plots
label <- c(
  "Trait aggression", 
  "Machiavellianism", 
  "Psychopathy", 
  "Narcissism", 
  "Perceived inability to gain status", 
  "Perceived status loss", 
  "Perceived marginalization", 
  "Authoritarian aggression", 
  "Populist attitudes"
)

color <- c(
  "#FDE725FF",  
  "#DCE319FF", 
  "#C7E020FF", 
  "#8FD744FF", 
  "#20A387FF", 
  "#2D706EFF", 
  "#31688DFF", 
  "#404788FF", 
  "#440154FF"
)

# Fit various OLS models (one per measured independent variable) with clustering by respondent ID (`RESPID`)

MB1.1_OLS_FIT <- lm.cluster(
  DV_UNCIVIL ~ POLATT_POPULISM + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES)

MB1.2_OLS_FIT <- lm.cluster(
  DV_UNCIVIL ~ POLATT_AUTHORITY + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES)

MB1.3_OLS_FIT <- lm.cluster(
  DV_UNCIVIL ~ SOCEXP_MARGINALZ + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES)

MB1.4_OLS_FIT <- lm.cluster(
  DV_UNCIVIL ~ SOCEXP_STATSLOSS + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES)

MB1.5_OLS_FIT <- lm.cluster(
  DV_UNCIVIL ~ SOCEXP_STATSGAIN + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES)

MB1.6_OLS_FIT <- lm.cluster(
  DV_UNCIVIL ~ PERSONL_NARRCISM + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES)

MB1.7_OLS_FIT <- lm.cluster(
  DV_UNCIVIL ~ PERSONL_PSYCHOPT + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES)

MB1.8_OLS_FIT <- lm.cluster(
  DV_UNCIVIL ~ PERSONL_MACCHIVL + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES)

MB1.9_OLS_FIT <- lm.cluster(
  DV_UNCIVIL ~ PERSONL_AGGRESSV + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES)

# Print a summary of the results for all models (Appendix E, Table E3)
export_summs(
  MB1.1_OLS_FIT$lm_res, MB1.2_OLS_FIT$lm_res, MB1.3_OLS_FIT$lm_res, 
  MB1.4_OLS_FIT$lm_res, MB1.5_OLS_FIT$lm_res, MB1.6_OLS_FIT$lm_res, 
  MB1.7_OLS_FIT$lm_res, MB1.8_OLS_FIT$lm_res, MB1.9_OLS_FIT$lm_res,
             coefs= c("Populist attitudes" = "POLATT_POPULISM",
                      "Authoritarian aggression" = "POLATT_AUTHORITY",
                      "Perceived marginalization" = "SOCEXP_MARGINALZ",
                      "Perceived status loss" = "SOCEXP_STATSLOSS",
                      "Perceived inability to gain status" = "SOCEXP_STATSGAIN",
                      "Narcissism" = "PERSONL_NARRCISM",
                      "Psychopathy" = "PERSONL_PSYCHOPT",
                      "Machiavellianism" = "PERSONL_MACCHIVL",
                      "Aggressiveness" = "PERSONL_AGGRESSV",
                      "Male" = "SEXMale", "Age" = "AGE", "Education" = "EDU_NUM",
                      "Political interest" = "as.numeric(POLINT)", 
                      "Party ID" = "as.numeric(PID_STR)", 
                      "Issue importance" = "as.numeric(ISSUE_IMPORTANCE)", 
                      "Outparty" = "ATT_PARTYDIROutparty", 
                      "Inparty" = "ATT_PARTYDIRInparty", "(Intercept)")
  )

# Create and display a forest plot for the model results (manuscript, Figure 3)
# This plot illustrates the estimates from all models for visual comparison
fun_forest(
  list(
    MB1.1_OLS_FIT$lm_res, MB1.2_OLS_FIT$lm_res, MB1.3_OLS_FIT$lm_res, 
    MB1.4_OLS_FIT$lm_res, MB1.5_OLS_FIT$lm_res, MB1.6_OLS_FIT$lm_res, 
    MB1.7_OLS_FIT$lm_res, MB1.8_OLS_FIT$lm_res, MB1.9_OLS_FIT$lm_res), 
  label,
  color,
  c("MB1.1", "MB1.2", "MB1.3", "MB1.4", "MB1.5", "MB1.6", "MB1.7", "MB1.8", "MB1.9"),
  "Estimates", 
  title2 = "")

## ============ Section 4.1: Robustness check 1; Incivility type (Models B2, B3, Appendix E, Figure E3) ============

# Code for robustness check 1: replication of MB1 per incivility type (Model B2: "insult"; Model B3: "ridicule") .
# Results are presented in appendix (Appendix E, Figure E3).

# Filter the dataset to include only insults and ridicule
DAT.ANALYSES.INS <- DAT.ANALYSES %>% filter(ATT_IMPOL == "Insult")
DAT.ANALYSES.RID <- DAT.ANALYSES %>% filter(ATT_IMPOL == "Ridicule")

# Fit OLS regression models for the "Insult" subset for each individual level predictors
# Model for Populist Attitudes
fit.polatt_popp_ins <- lm.cluster(
  DV_UNCIVIL ~ POLATT_POPULISM + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.INS
)

# Model for Authoritarian Aggression
fit.polatt_auth_ins <- lm.cluster(
  DV_UNCIVIL ~ POLATT_AUTHORITY + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.INS
)

# Model for Perceived Marginalization
fit.socexp_marg_ins <- lm.cluster(
  DV_UNCIVIL ~ SOCEXP_MARGINALZ + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.INS
)

# Model for Perceived Status Loss
fit.socexp_loss_ins <- lm.cluster(
  DV_UNCIVIL ~ SOCEXP_STATSLOSS + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.INS
)

# Model for Perceived Inability to Gain Status
fit.socexp_gain_ins <- lm.cluster(
  DV_UNCIVIL ~ SOCEXP_STATSGAIN + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.INS
)

# Model for Narcissism
fit.person_narc_ins <- lm.cluster(
  DV_UNCIVIL ~ PERSONL_NARRCISM + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.INS
)

# Model for Psychopathy
fit.person_pych_ins <- lm.cluster(
  DV_UNCIVIL ~ PERSONL_PSYCHOPT + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.INS
)

# Model for Machiavellianism
fit.person_mach_ins <- lm.cluster(
  DV_UNCIVIL ~ PERSONL_MACCHIVL + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.INS
)

# Model for Aggressiveness
fit.person_aggr_ins <- lm.cluster(
  DV_UNCIVIL ~ PERSONL_AGGRESSV + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.INS
)

# Fit OLS regression models for the "Ridicule" subset for each individual level predictors

# Model for Populist Attitudes
fit.polatt_popp_rid <- lm.cluster(
  DV_UNCIVIL ~ POLATT_POPULISM + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.RID
)

# Model for Authoritarian Aggression
fit.polatt_auth_rid <- lm.cluster(
  DV_UNCIVIL ~ POLATT_AUTHORITY + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.RID
)

# Model for Perceived Marginalization
fit.socexp_marg_rid <- lm.cluster(
  DV_UNCIVIL ~ SOCEXP_MARGINALZ + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.RID
)

# Model for Perceived Status Loss
fit.socexp_loss_rid <- lm.cluster(
  DV_UNCIVIL ~ SOCEXP_STATSLOSS + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.RID
)

# Model for Perceived Inability to Gain Status
fit.socexp_gain_rid <- lm.cluster(
  DV_UNCIVIL ~ SOCEXP_STATSGAIN + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.RID
)

# Model for Narcissism
fit.person_narc_rid <- lm.cluster(
  DV_UNCIVIL ~ PERSONL_NARRCISM + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.RID
)

# Model for Psychopathy
fit.person_pych_rid <- lm.cluster(
  DV_UNCIVIL ~ PERSONL_PSYCHOPT + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.RID
)

# Model for Machiavellianism
fit.person_mach_rid <- lm.cluster(
  DV_UNCIVIL ~ PERSONL_MACCHIVL + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.RID
)

# Model for Aggressiveness
fit.person_aggr_rid <- lm.cluster(
  DV_UNCIVIL ~ PERSONL_AGGRESSV + AGE + SEX + EDU_NUM + 
    as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.ANALYSES.RID
)


# Generate forest plots for the "Insult" models
MB2_OLS_PLOT <- fun_forest(
  list(
    fit.polatt_popp_ins$lm_res, fit.polatt_auth_ins$lm_res, fit.socexp_marg_ins$lm_res, 
    fit.socexp_loss_ins$lm_res, fit.socexp_gain_ins$lm_res, fit.person_narc_ins$lm_res,
    fit.person_pych_ins$lm_res, fit.person_mach_ins$lm_res, fit.person_aggr_ins$lm_res
    ), 
  label,
  color,
  c("MB2.1", "MB2.2", "MB2.3", "MB2.4", "MB2.5", "MB2.6", "MB2.7", "MB2.8", "MB2.9"),
  title ="Estimates",
  title2 = "Insulting language") + 
  geom_point(shape = 1, size = 3)
                           

# Generate forest plots for the "Ridicule" models
MB3_OLS_PLOT <-fun_forest(
  list(fit.polatt_popp_rid$lm_res, fit.polatt_auth_rid$lm_res, fit.socexp_marg_rid$lm_res, 
       fit.socexp_loss_rid$lm_res, fit.socexp_gain_rid$lm_res, fit.person_narc_rid$lm_res, 
       fit.person_pych_rid$lm_res, fit.person_mach_rid$lm_res, fit.person_aggr_rid$lm_res),
  label,
  color,
  c("MB3.1", "MB3.2", "MB3.3", "MB3.4", "MB3.5", "MB3.6", "MB3.7", "MB3.8", "MB3.9"),
  title ="Estimates",
  title2 = "Ridicule") + 
  geom_point(shape = 2, size = 3)
                          

# Combine the forest plots for "Insult" and "Ridicule" into a single plot layout (Appendix E, Figure E3)
egg::ggarrange(
  MB2_OLS_PLOT + theme(
    legend.position = "none", legend.justification = "left"
    ), 
  MB3_OLS_PLOT + theme(
    legend.position = "none", legend.justification = "left", axis.text.y = element_blank()
    ), 
  nrow = 1)

## ============ Section 4.2: Robustness check 2; Impoliteness & inappropriateness as DVs (B4, B5, B6, B7, Appendix E, Figure E4) ============

# Code for robustness check 2: replication of MB1 per incivility perceptions indicators
# Model B4: "impolite"; Model B5: "wrong"; Model B6: "unacceptable"; Model B7: "not normal".
# Results are presented in appendix (Appendix E, Figure E4).

# Fit linear models with clustered standard errors for each individual level predictors and outcome variables
# 1. Fit models for the "Impolite" dependent variable
fit.polatt_popp_imp <- lm.cluster(
  DV_IMPOLITE ~ POLATT_POPULISM + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.polatt_auth_imp <- lm.cluster(
  DV_IMPOLITE ~ POLATT_AUTHORITY + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.socexp_marg_imp <- lm.cluster(
  DV_IMPOLITE ~ SOCEXP_MARGINALZ + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.socexp_loss_imp <- lm.cluster(
  DV_IMPOLITE ~ SOCEXP_STATSLOSS + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.socexp_gain_imp <- lm.cluster(
  DV_IMPOLITE ~ SOCEXP_STATSGAIN + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_narc_imp <- lm.cluster(
  DV_IMPOLITE ~ PERSONL_NARRCISM + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_pych_imp <- lm.cluster(
  DV_IMPOLITE ~ PERSONL_PSYCHOPT + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_mach_imp <- lm.cluster(
  DV_IMPOLITE ~ PERSONL_MACCHIVL + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_aggr_imp <- lm.cluster(
  DV_IMPOLITE ~ PERSONL_AGGRESSV + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

# 2. Fit models for the "Wrong" dependent variable
fit.polatt_popp_wrong <- lm.cluster(
  DV_WRONG ~ POLATT_POPULISM + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.polatt_auth_wrong <- lm.cluster(
  DV_WRONG ~ POLATT_AUTHORITY + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.socexp_marg_wrong <- lm.cluster(
  DV_WRONG ~ SOCEXP_MARGINALZ + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.socexp_loss_wrong <- lm.cluster(
  DV_WRONG ~ SOCEXP_STATSLOSS + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.socexp_gain_wrong <- lm.cluster(
  DV_WRONG ~ SOCEXP_STATSGAIN + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_narc_wrong <- lm.cluster(
  DV_WRONG ~ PERSONL_NARRCISM + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_pych_wrong <- lm.cluster(
  DV_WRONG ~ PERSONL_PSYCHOPT + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_mach_wrong <- lm.cluster(
  DV_WRONG ~ PERSONL_MACCHIVL + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_aggr_wrong <- lm.cluster(
  DV_WRONG ~ PERSONL_AGGRESSV + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

# 3. Fit models for the "Unacceptable" dependent variable
fit.polatt_popp_unacp <- lm.cluster(
  DV_UNACCEPT ~ POLATT_POPULISM + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.polatt_auth_unacp <- lm.cluster(
  DV_UNACCEPT ~ POLATT_AUTHORITY + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.socexp_marg_unacp <- lm.cluster(
  DV_UNACCEPT ~ SOCEXP_MARGINALZ + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.socexp_loss_unacp <- lm.cluster(
  DV_UNACCEPT ~ SOCEXP_STATSLOSS + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.socexp_gain_unacp <- lm.cluster(
  DV_UNACCEPT ~ SOCEXP_STATSGAIN + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_narc_unacp <- lm.cluster(
  DV_UNACCEPT ~ PERSONL_NARRCISM + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_pych_unacp <- lm.cluster(
  DV_UNACCEPT ~ PERSONL_PSYCHOPT + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_mach_unacp <- lm.cluster(
  DV_UNACCEPT ~ PERSONL_MACCHIVL + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_aggr_unacp <- lm.cluster(
  DV_UNACCEPT ~ PERSONL_AGGRESSV + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

# 4. Fit models for the "Not normal" dependent variable
fit.polatt_popp_norml <- lm.cluster(
  DV_NOTNORMAL ~ POLATT_POPULISM + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.polatt_auth_norml <- lm.cluster(
  DV_NOTNORMAL ~ POLATT_AUTHORITY + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.socexp_marg_norml <- lm.cluster(
  DV_NOTNORMAL ~ SOCEXP_MARGINALZ + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.socexp_loss_norml <- lm.cluster(
  DV_NOTNORMAL ~ SOCEXP_STATSLOSS + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.socexp_gain_norml <- lm.cluster(
  DV_NOTNORMAL ~ SOCEXP_STATSGAIN + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_narc_norml <- lm.cluster(
  DV_NOTNORMAL ~ PERSONL_NARRCISM + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_pych_norml <- lm.cluster(
  DV_NOTNORMAL ~ PERSONL_PSYCHOPT + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_mach_norml <- lm.cluster(
  DV_NOTNORMAL ~ PERSONL_MACCHIVL + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

fit.person_aggr_norml <- lm.cluster(
  DV_NOTNORMAL ~ PERSONL_AGGRESSV + AGE + SEX + EDU_NUM + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID", 
  data = DAT.ANALYSES
)

# Create a forest plot for each set of model

# Creates a forest plot for models with DV_IMPOLITE
MB4_OLS_PLOT <- fun_forest(
  list(
    fit.polatt_popp_imp$lm_res, 
    fit.polatt_auth_imp$lm_res, 
    fit.socexp_marg_imp$lm_res, 
    fit.socexp_loss_imp$lm_res, 
    fit.socexp_gain_imp$lm_res, 
    fit.person_narc_imp$lm_res, 
    fit.person_pych_imp$lm_res, 
    fit.person_mach_imp$lm_res, 
    fit.person_aggr_imp$lm_res
  ), 
  label,
  color,
  c("MB4.1", "MB4.2", "MB4.3", "MB4.4", "MB4.5", "MB4.6", "MB4.7", "MB4.8", "MB4.9"),
  title = "Estimates",
  title2 = "DV: Impolite"
) + geom_point(shape = 3, size = 3)

# Creates a forest plot for models with DV_WRONG
MB5_OLS_PLOT <- fun_forest(
  list(
    fit.polatt_popp_wrong$lm_res, 
    fit.polatt_auth_wrong$lm_res, 
    fit.socexp_marg_wrong$lm_res, 
    fit.socexp_loss_wrong$lm_res, 
    fit.socexp_gain_wrong$lm_res, 
    fit.person_narc_wrong$lm_res, 
    fit.person_pych_wrong$lm_res, 
    fit.person_mach_wrong$lm_res, 
    fit.person_aggr_wrong$lm_res
  ), 
  label,
  color,
  c("MB5.1", "MB5.2", "MB5.3", "MB5.4", "MB5.5", "MB5.6", "MB5.7", "MB5.8", "MB5.9"),
  title = "Estimates",
  title2 = "DV: Wrong"
) + geom_point(shape = 4, size = 3)

# Creates a forest plot for models with DV_UNACCEPT
MB6_OLS_PLOT <- fun_forest(
  list(
    fit.polatt_popp_unacp$lm_res, 
    fit.polatt_auth_unacp$lm_res, 
    fit.socexp_marg_unacp$lm_res, 
    fit.socexp_loss_unacp$lm_res, 
    fit.socexp_gain_unacp$lm_res, 
    fit.person_narc_unacp$lm_res, 
    fit.person_pych_unacp$lm_res, 
    fit.person_mach_unacp$lm_res, 
    fit.person_aggr_unacp$lm_res
  ), 
  label,
  color,
  c("MB6.1", "MB6.2", "MB6.3", "MB6.4", "MB6.5", "MB6.6", "MB6.7", "MB6.8", "MB6.9"),
  title = "Estimates",
  title2 = "DV: Unacceptable"
) + geom_point(shape = 5, size = 3)

# Creates a forest plot for models with DV_NOTNORMAL
MB7_OLS_PLOT <- fun_forest(
  list(
    fit.polatt_popp_norml$lm_res, 
    fit.polatt_auth_norml$lm_res, 
    fit.socexp_marg_norml$lm_res, 
    fit.socexp_loss_norml$lm_res, 
    fit.socexp_gain_norml$lm_res, 
    fit.person_narc_norml$lm_res, 
    fit.person_pych_norml$lm_res, 
    fit.person_mach_norml$lm_res, 
    fit.person_aggr_norml$lm_res
  ), 
  label,
  color,
  c("MB7.1", "MB7.2", "MB7.3", "MB7.4", "MB7.5", "MB7.6", "MB7.7", "MB7.8", "MB7.9"),
  title = "Estimates",
  title2 = "DV: Not normal"
) + geom_point(shape = 6, size = 3)

# Combine multiple forest plots into a single layout (Appendix E, Figure E4)
egg::ggarrange(
  MB4_OLS_PLOT + theme(legend.position = "none"), 
  MB5_OLS_PLOT + theme(legend.position = "none", axis.text.y = element_blank()), 
  MB6_OLS_PLOT + theme(legend.position = "none", axis.text.y = element_blank()),
  MB7_OLS_PLOT + theme(legend.position = "none", axis.text.y = element_blank()), 
  nrow = 1
)

# Combine multiple forest plots into a single layout (manuscript, Figure 4)
egg::ggarrange(
  MB7_OLS_PLOT + theme(legend.position = "none"), 
  MB2_OLS_PLOT + theme(legend.position = "none", axis.text.y = element_blank()), 
  MB3_OLS_PLOT + theme(legend.position = "none", axis.text.y = element_blank()), 
  nrow = 1
)

## ============ Section 4.3: Sensitivity Analysis (Appendix E) ============

# Result from the sensitivity analysis for Model B1, as reported in Appendix E

# Function to perform sensitivity analysis using the sensemakr package
perform_sensitivity_analysis <- function(model, treatment, benchmark_covariates, kd) {
  sensemakr(
    model = model,
    treatment = treatment,
    benchmark_covariates = benchmark_covariates,
    kd = kd
  )
}

# Sensitivity analysis for Model B1 as reported in Appendix E
perform_sensitivity_analysis(
  model = MB1.1_OLS_FIT$lm_res,
  treatment = "POLATT_POPULISM",
  benchmark_covariates = "ATT_PARTYDIROutparty",
  kd = 1:3
)

perform_sensitivity_analysis(
  model = MB1.2_OLS_FIT$lm_res,
  treatment = "POLATT_AUTHORITY",
  benchmark_covariates = "ATT_PARTYDIROutparty",
  kd = 1:3
)

perform_sensitivity_analysis(
  model = MB1.3_OLS_FIT$lm_res,
  treatment = "SOCEXP_MARGINALZ",
  benchmark_covariates = "ATT_PARTYDIROutparty",
  kd = 1:3
)

perform_sensitivity_analysis(
  model = MB1.6_OLS_FIT$lm_res,
  treatment = "PERSONL_NARRCISM",
  benchmark_covariates = "ATT_PARTYDIROutparty",
  kd = 1:3
)

perform_sensitivity_analysis(
  model = MB1.7_OLS_FIT$lm_res,
  treatment = "PERSONL_PSYCHOPT",
  benchmark_covariates = "ATT_PARTYDIROutparty",
  kd = 1:3
)

perform_sensitivity_analysis(
  model = MB1.8_OLS_FIT$lm_res,
  treatment = "PERSONL_MACCHIVL",
  benchmark_covariates = "ATT_PARTYDIROutparty",
  kd = 1:3
)

perform_sensitivity_analysis(
  model = MB1.9_OLS_FIT$lm_res,
  treatment = "PERSONL_AGGRESSV",
  benchmark_covariates = "ATT_PARTYDIROutparty",
  kd = 1:3
)

# ============ Section 5: Exploratory Analysis; Mediation models (Model C1) ============

# This section conducts the exploratory analyses on mediation effect of incivility perceptions for Model C1.
# The results are presented in the manuscript (pp. 22-23, Figure 5 & 6) and in the appendix (Appendix E, Table E4 & E5).
# Note: The results of the mediation models vary slightly in each run due to Monte Carlo sampling errors.
# This occurs because the mediate() function uses bootstrapping and Monte Carlo simulations to estimate indirect effects.
# To ensure reproducibility of results, we set a random seed (set.seed(1234)) before each iteration.

# Disable scientific notation for better readability
options(scipen = 999)


# Select relevant variables and remove missing values
DAT.MED <- DAT.ANALYSES %>% 
  dplyr::select(
    DV_FEELING, DV_UNCIVIL, ATT_GENDR, ATT_TRAIT, ATT_POPUL, ATT_ATTAK, 
    ATT_ISSUE, ATT_SETTG, SEX, POLINT, PID_STR, ISSUE_IMPORTANCE, 
    ATT_PARTYDIR, RESPID
  ) %>% 
  drop_na()

# Step 1) Fit linear models for the mediator variable (DV_UNCIVIL) per selected attribute

# fit for attribute: personality/leadership traits
fit.med.trait <- lm.cluster(
  DV_UNCIVIL ~ ATT_TRAIT + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.MED
)

# fit for attribute: incivility setting
fit.med.settg <- lm.cluster(
  DV_UNCIVIL ~ ATT_SETTG + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.MED
)

# fit for attribute: partisan congruency
fit.med.party <- lm.cluster(
  DV_UNCIVIL ~ ATT_PARTYDIR + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE), 
  cluster = "RESPID", 
  data = DAT.MED
)

# Step 2) Fit linear models for the outcome variable (DV_FEELING) clustered per selected attribute

# fit for attribute: personality/leadership traits
fit.out.trait <- lm.cluster(
  DV_FEELING ~ ATT_TRAIT + DV_UNCIVIL + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.MED
)

# fit for attribute: incivility setting
fit.out.settg <- lm.cluster(
  DV_FEELING ~ ATT_SETTG + DV_UNCIVIL + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.MED
)

# fit for attribute: partisan congruency
fit.out.party <- lm.cluster(
  DV_FEELING ~ ATT_PARTYDIR + DV_UNCIVIL + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE), 
  cluster = "RESPID", 
  data = DAT.MED
)

# Step 3) Perform mediation analysis using the fitted models per selected attribute

# fit for attribute: personality/leadership traits
set.seed(1234)
medmodel.trait <- mediation::mediate(
  fit.med.trait$lm_res, 
  fit.out.trait$lm_res, 
  treat = "ATT_TRAIT", 
  mediator = "DV_UNCIVIL", 
  robustSE = TRUE, 
  sims = 1000
)

# fit for attribute: incivility settings
set.seed(1234)
medmodel.settg <- mediation::mediate(
  fit.med.settg$lm_res, 
  fit.out.settg$lm_res, 
  treat = "ATT_SETTG", 
  mediator = "DV_UNCIVIL", 
  robustSE = TRUE, 
  sims = 1000
)

# fit for attribute: partisanship (outparty vs unknown)
set.seed(1234)
medmodel.party1 <- mediation::mediate(
  fit.med.party$lm_res, 
  fit.out.party$lm_res, 
  treat = "ATT_PARTYDIR", 
  mediator = "DV_UNCIVIL", 
  control.value = "Unknown", 
  treat.value = "Outparty", 
  robustSE = TRUE, 
  sims = 1000
)

# fit for attribute: partisanship (inparty vs unknown)
set.seed(1234)
medmodel.party2 <- mediation::mediate(
  fit.med.party$lm_res, 
  fit.out.party$lm_res, 
  treat = "ATT_PARTYDIR", 
  mediator = "DV_UNCIVIL", 
  control.value = "Unknown", 
  treat.value = "Inparty", 
  robustSE = TRUE, 
  sims = 1000
)

# fit for attribute: partisanship (inparty vs outparty)
set.seed(1234)
medmodel.party3 <- mediation::mediate(
  fit.med.party$lm_res, 
  fit.out.party$lm_res, 
  treat = "ATT_PARTYDIR", 
  mediator = "DV_UNCIVIL", 
  control.value = "Outparty", 
  treat.value = "Inparty", 
  robustSE = TRUE, 
  sims = 1000
)

# Step 4) Summarize results from the mediation models

summary(medmodel.trait) # Model C1.1, Appendix E (Table E4)
summary(medmodel.settg) # Model C1.2, Appendix E (Table E4)
summary(medmodel.party1) # Model C1.3, Appendix E (Table E5)
summary(medmodel.party2) # Model C1.4, Appendix E (Table E5)
summary(medmodel.party3) # Model C1.5, Appendix E (Table E5)

# Step 5) Plot results for mediation models
plot(medmodel.trait, main = "Leadership Traits", xlim = c(-10, 2))  # Plot for personality/leadership traits (Figure 5)
plot(medmodel.settg, main = "Incivility Setting", xlim = c(-10, 2))  # Plot for incivility setting (Figure 5)
plot(medmodel.party1, main = "Outparty (Ref: Unknown)", xlim = c(-20, 20))  # Plot for partisanship (Figure 6)
plot(medmodel.party2, main = "Inparty (Ref: Unknown)", xlim = c(-20, 20))  # Plot for partisanship (Figure 6)
plot(medmodel.party3, main = "Inparty (Ref: Outparty)", xlim = c(-20, 20))  # Plot for partisanship (Figure 6)

## ============ Section 5.1: Robustness check 1; Incivility type (Models C2.1 - C2.4, C3.1 - C3.4) ============

# Code for Robustness Check 1: Replication of MC1 by Incivility Type
# Models: C2 ("insult") and C3 ("ridicule")
# Results are for Appendix E, Table E6 and Table E7

#### ============ Section 5.1.1 Mediation for "Insult" models (MC2.1, MC2.2, Table E6) ============

# Filter dataset to include only "insults" 
DAT.ANALYSES.INS <- DAT.ANALYSES %>% 
  filter(ATT_IMPOL == "Insult")

# Select relevant variables for mediation analysis and remove missing values
DAT.MED.INS <- DAT.ANALYSES.INS %>% 
  dplyr::select(
    DV_FEELING, DV_UNCIVIL, ATT_GENDR, ATT_TRAIT, ATT_POPUL, ATT_ATTAK, 
    ATT_ISSUE, ATT_SETTG, SEX, POLINT, PID_STR, ISSUE_IMPORTANCE, 
    ATT_PARTYDIR, RESPID
    ) %>% 
  drop_na()

# Step 1) Fit linear models for the mediator variable (DV_UNCIVIL), clustered by respondent ID (RESPID)

# fit for attribute: personality/leadership traits
fit.med.trait.ins <- lm.cluster(
  DV_UNCIVIL ~ ATT_TRAIT + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.MED.INS)

# fit for attribute: incivility setting
fit.med.settg.ins <- lm.cluster(
  DV_UNCIVIL ~ ATT_SETTG + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.MED.INS)

# Step 2) Fit linear models for the outcome variable (DV_FEELING), clustered by respondent ID (RESPID)

# fit for attribute: personality/leadership traits
fit.out.trait.ins <- lm.cluster(
  DV_FEELING ~ ATT_TRAIT + DV_UNCIVIL + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.MED.INS)

# fit for attribute: incivility setting
fit.out.settg.ins <- lm.cluster(
  DV_FEELING ~ ATT_SETTG + DV_UNCIVIL + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.MED.INS)

# Step 3) Mediation analysis testing the role of DV_UNCIVIL as a mediator

# fit for attribute: personality/leadership traits
set.seed(1234)
medmodel.trait.ins <- mediation::mediate(
  fit.med.trait.ins$lm_res, 
  fit.out.trait.ins$lm_res, 
  treat = "ATT_TRAIT", 
  mediator = "DV_UNCIVIL", 
  robustSE = TRUE, 
  sims = 1000)

# fir for attribute: incivility setting
set.seed(1234)
medmodel.settg.ins <- mediation::mediate(
  fit.med.settg.ins$lm_res, 
  fit.out.settg.ins$lm_res, 
  treat = "ATT_SETTG", 
  mediator = "DV_UNCIVIL", 
  robustSE = TRUE, 
  sims = 1000)

# Summarize mediation results for "insult" (Appendix E, Table E6)
summary(medmodel.trait.ins) # attribute: personality/leadership traits (Model C2.1)
summary(medmodel.settg.ins) # attribute: incivility setting (Model C2.2)

#### ============ Section 5.1.2 Mediation for "Ridicule" models (MC3.1, MC3.2, Table E7) ============

# Filter dataset to include only "ridicule"
DAT.ANALYSES.RID <- DAT.ANALYSES %>% 
  filter(ATT_IMPOL == "Ridicule")

# Select relevant variables for mediation analysis and remove missing values
DAT.MED.RID <- DAT.ANALYSES.RID %>% 
  dplyr::select(
    DV_FEELING, DV_UNCIVIL, ATT_GENDR, ATT_TRAIT, ATT_POPUL, ATT_ATTAK, 
    ATT_ISSUE, ATT_SETTG, SEX, POLINT, PID_STR, ISSUE_IMPORTANCE, 
    ATT_PARTYDIR, RESPID
    ) %>% 
  drop_na()

# Step 1) Fit linear models for the mediator variable (DV_UNCIVIL) per selected attribute

# fit for attribute: personality/leadership traits
fit.med.trait.rid <- lm.cluster(
  DV_UNCIVIL ~ ATT_TRAIT + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.MED.RID)

# attribute: incivility setting
fit.med.settg.rid <- lm.cluster(
  DV_UNCIVIL ~ ATT_SETTG + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.MED.RID)

# Step 2) Fit linear models for the outcome variable (DV_FEELING) per selected attribute

# fit for attribute: personality/leadership traits
fit.out.trait.rid <- lm.cluster(
  DV_FEELING ~ ATT_TRAIT + DV_UNCIVIL + SEX + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.MED.RID)

# fit for attribute: incivility setting
fit.out.settg.rid <- lm.cluster(
  DV_FEELING ~ ATT_SETTG + DV_UNCIVIL + SEX + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data = DAT.MED.RID)

# Step 3) Mediation analysis testing the role of DV_UNCIVIL as a mediator per selected attribute

# attribute: personality/leadership traits
set.seed(1234)
medmodel.trait.rid <- mediation::mediate(
  fit.med.trait.rid$lm_res, 
  fit.out.trait.rid$lm_res, 
  treat = "ATT_TRAIT", 
  mediator = "DV_UNCIVIL", 
  robustSE = TRUE, 
  sims = 1000)

# attribute: incivility setting ()
set.seed(1234)
medmodel.settg.rid <- mediation::mediate(
  fit.med.settg.rid$lm_res, 
  fit.out.settg.rid$lm_res, 
  treat = "ATT_SETTG", mediator = "DV_UNCIVIL", 
  robustSE = TRUE, 
  sims = 1000)

# Step 4) Summarize mediation results for "ridicule" (Appendix E, Table E7)
summary(medmodel.trait.rid) # attribute: personality/leadership traits (Model C3.1)
summary(medmodel.settg.rid) # attribute: incivility setting (Model C3.2)


## ============ Section 5.2: Robustness check 2; Impoliteness & inappropriateness as DVs (Model C4, C5, C6, C7, C8) ============

# Code for robustness check 2: replication of MC1 per incivility perceptions indicators
# Model C4: "impolite"; Model C5: "wrong"; Model C6: "unacceptable"; Model C7: "not normal".
# Results for appendix: Appendix E, Figure E8 - Figure E11.

#### ============ Section 5.2.1: Mediation models with DV_IMPOLITE (Model C4, Appendix E, Table E8) ============

options(scipen = 999)

# Select relevant variables for mediation analysis and remove missing values
DAT.MED.IMP <- DAT.ANALYSES %>% 
  dplyr::select(
    DV_FEELING, DV_IMPOLITE, ATT_GENDR, ATT_TRAIT, ATT_POPUL, ATT_ATTAK, 
    ATT_ISSUE, ATT_SETTG, SEX, POLINT, PID_STR, ISSUE_IMPORTANCE, 
    ATT_PARTYDIR, RESPID
    ) %>% 
  drop_na()

# Step 1) Fit linear models for the mediator variable (DV_IMPOLITE) per selected attributes

# fit for attribute: personality/leadership traits
fit.med.trait.imp <- lm.cluster(
  DV_IMPOLITE ~ ATT_TRAIT + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data =  DAT.MED.IMP)

# fit for attribute: incivility settings
fit.med.settg.imp <- lm.cluster(
  DV_IMPOLITE ~ ATT_SETTG + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data =  DAT.MED.IMP)

# Step 2) Fit linear models for the outcome variable (DV_FEELING) per selected attributes

# fit for attribute: personality/leadership traits
fit.out.trait.imp <- lm.cluster(
  DV_FEELING ~ ATT_TRAIT + DV_IMPOLITE + SEX + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data =  DAT.MED.IMP)

# fit for attribute: incivility settings
fit.out.settg.imp <- lm.cluster(
  DV_FEELING ~ ATT_SETTG + DV_IMPOLITE + SEX + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data =  DAT.MED.IMP)

# Step 3) Mediation analysis testing the role of DV_IMPOLITE as a mediator per selected attribute

# fit for attribute: personality/leadership traits
set.seed(1234)
medmodel.trait.imp <- mediation::mediate(
  fit.med.trait.imp$lm_res, 
  fit.out.trait.imp$lm_res, 
  treat = "ATT_TRAIT", 
  mediator = "DV_IMPOLITE", 
  robustSE = TRUE, 
  sims = 1000)

# fit for attribute: incivility settings
set.seed(1234)
medmodel.settg.imp <- mediation::mediate(
  fit.med.settg.imp$lm_res, 
  fit.out.settg.imp$lm_res,
  treat = "ATT_SETTG", 
  mediator = "DV_IMPOLITE", 
  robustSE = TRUE, 
  sims = 1000)

# Step 4) Summarize mediation results (Appendix E, Table E8)
summary(medmodel.trait.imp) # attribute: personality/leadership traits (Model C4.1)
summary(medmodel.settg.imp) # attribute: incivility setting (Model C4.2)


#### ============ Section 5.2.2: Mediation models with DV_WRONG (Model C5, Appendix E, Table E9) ============

options(scipen = 999)

# Select relevant variables for mediation analysis and remove missing values
DAT.MED.WRG <- DAT.ANALYSES %>% 
  dplyr::select(
    DV_FEELING, DV_WRONG, ATT_GENDR, ATT_TRAIT, ATT_POPUL, ATT_ATTAK, 
    ATT_ISSUE, ATT_SETTG, SEX, POLINT, PID_STR, ISSUE_IMPORTANCE, 
    ATT_PARTYDIR, RESPID
    ) %>% 
  drop_na()

# Step 1) Fit linear models for the mediator variable (DV_WRONG) per selected attribute

# fit for attribute: personality/leadership traits
fit.med.trait.wrg <- lm.cluster(
  DV_WRONG ~ ATT_TRAIT + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID",
  data = DAT.MED.WRG
)

# fit for attribute: incivility settings
fit.med.settg.wrg <- lm.cluster(
  DV_WRONG ~ ATT_SETTG + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID",
  data = DAT.MED.WRG
  )

# Step 2) Fit linear models for the outcome variable (DV_FEELING) per selected attribute

# fit for attribute: personality/leadership traits
fit.out.trait.wrg <- lm.cluster(
  DV_FEELING ~ ATT_TRAIT + DV_WRONG + SEX + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID",
  data = DAT.MED.WRG
)

# fit for attribute: incivility settings
fit.out.settg.wrg <- lm.cluster(
  DV_FEELING ~ ATT_SETTG + DV_WRONG + SEX + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR,
  cluster = "RESPID",
  data = DAT.MED.WRG
)

# Step 3) Mediation analysis testing the role of DV_WRONG as a mediator per selected attribute

# fit for attribute: personality/leadership traits
set.seed(1234)
medmodel.trait.wrg <- mediation::mediate(
  fit.med.trait.wrg$lm_res,
  fit.out.trait.wrg$lm_res,
  treat = "ATT_TRAIT",
  mediator = "DV_WRONG",
  robustSE = TRUE,
  sims = 1000
)

# fir for attribute: incivility settings
set.seed(1234)
medmodel.settg.wrg <- mediation::mediate(
  fit.med.settg.wrg$lm_res,
  fit.out.settg.wrg$lm_res,
  treat = "ATT_SETTG",
  mediator = "DV_WRONG",
  robustSE = TRUE,
  sims = 1000
)

# Step 4) Summarize mediation results (Appendix E, Table E9)
summary(medmodel.trait.wrg) # attribute: personality/leadership traits (Model C5.1)
summary(medmodel.settg.wrg) # attribute: incivility setting (Model C5.2)


#### ============ Section 5.2.3: Mediation models with DV_UNACCEPTABLE (Model C6, Appendix E, Table E10) ============

options(scipen = 999)

# Select relevant variables for mediation analysis and remove missing values
DAT.MED.INP <- DAT.ANALYSES %>% 
  dplyr::select(
    DV_FEELING, DV_UNACCEPT, ATT_GENDR, ATT_TRAIT, ATT_POPUL, ATT_ATTAK, 
    ATT_ISSUE, ATT_SETTG, SEX, POLINT, PID_STR, ISSUE_IMPORTANCE, 
    ATT_PARTYDIR, RESPID) %>% 
  drop_na()

# Step 1) Fit linear models for the mediator variable (DV_UNACCEPT) per selected attribute

# fit for attribute: personality/leadership traits
fit.med.trait.inp <- lm.cluster(
  DV_UNACCEPT ~ ATT_TRAIT + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data =  DAT.MED.INP)

# fit for attribute: incivility setting
fit.med.settg.inp <- lm.cluster(
  DV_UNACCEPT ~ ATT_SETTG + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data =  DAT.MED.INP)


# Step 2) Fit linear models for the outcome variable (DV_FEELING) per selected attributes

# fit for attribute: personality/leadership traits
fit.out.trait.inp <- lm.cluster(
  DV_FEELING ~ ATT_TRAIT + DV_UNACCEPT + SEX + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data =  DAT.MED.INP)

# fit for attribute: incivility setting
fit.out.settg.inp <- lm.cluster(
  DV_FEELING ~ ATT_SETTG + DV_UNACCEPT + SEX + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data =  DAT.MED.INP)

# Step 3) Mediation analysis testing the role of DV_WRONG as a mediator per selected attribute

# fit for attribute: personality/leadership traits
set.seed(1234)
medmodel.trait.inp <- mediation::mediate(
  fit.med.trait.inp$lm_res, 
  fit.out.trait.inp$lm_res, 
  treat = "ATT_TRAIT", 
  mediator = "DV_UNACCEPT", 
  robustSE = TRUE, 
  sims = 1000)

# fit for attribute: incivility setting
set.seed(1234)
medmodel.settg.inp <- mediation::mediate(
  fit.med.settg.inp$lm_res, 
  fit.out.settg.inp$lm_res, 
  treat = "ATT_SETTG", 
  mediator = "DV_UNACCEPT", 
  robustSE = TRUE, 
  sims = 1000)

# Step 4) Summarize mediation results (Appendix E, Table E10)
summary(medmodel.trait.inp) # attribute: personality/leadership traits (Model C6.1)
summary(medmodel.settg.inp) # attribute: incivility setting (Model C6.2)

#### ============ Section 5.2.4: Mediation models with DV_NOTNORMAL (Model C7, Appendix E, Table E11) ============

options(scipen = 999)

# Select relevant variables for mediation analysis and remove missing values
DAT.MED.NNR <- DAT.ANALYSES %>% 
  dplyr::select(
    DV_FEELING, DV_NOTNORMAL, ATT_GENDR, ATT_TRAIT, ATT_POPUL, ATT_ATTAK, 
    ATT_ISSUE, ATT_SETTG, SEX, POLINT, PID_STR, ISSUE_IMPORTANCE, 
    ATT_PARTYDIR, RESPID
    ) %>% 
  drop_na()

# Step 1) Fit linear models for the mediator variable (DV_NOTNORMAL) per selected attribute

# fit for attribute: personality/leadership traits
fit.med.trait.nnr <- lm.cluster(
  DV_NOTNORMAL ~ ATT_TRAIT + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data =  DAT.MED.NNR)

# fit for attribute: incivility setting
fit.med.settg.nnr <- lm.cluster(
  DV_NOTNORMAL ~ ATT_SETTG + SEX + as.numeric(POLINT) + as.numeric(PID_STR) + 
    as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data =  DAT.MED.NNR)

# Step 2) Fit linear models for the outcome variable (DV_FEELING) per selected attributes

# fit for attribute: personality/leadership traits
fit.out.trait.nnr <- lm.cluster(
  DV_FEELING ~ ATT_TRAIT + DV_NOTNORMAL + SEX + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data =  DAT.MED.NNR)

# fit for attribute: incivility setting
fit.out.settg.nnr <- lm.cluster(
  DV_FEELING ~ ATT_SETTG + DV_NOTNORMAL + SEX + as.numeric(POLINT) + 
    as.numeric(PID_STR) + as.numeric(ISSUE_IMPORTANCE) + ATT_PARTYDIR, 
  cluster = "RESPID", 
  data =  DAT.MED.NNR)

# Step 3) Mediation analysis testing the role of DV_WRONG as a mediator per selected attribute

# fit for attribute: personality/leadership traits
set.seed(1234)
medmodel.trait.nnr <- mediation::mediate(
  fit.med.trait.nnr$lm_res, 
  fit.out.trait.nnr$lm_res, 
  treat = "ATT_TRAIT", 
  mediator = "DV_NOTNORMAL", 
  robustSE = TRUE, 
  sims = 1000)

# fit for attribute: incivility setting
set.seed(1234)
medmodel.settg.nnr <- mediation::mediate(
  fit.med.settg.nnr$lm_res, 
  fit.out.settg.nnr$lm_res, 
  treat = "ATT_SETTG", 
  mediator = "DV_NOTNORMAL", 
  robustSE = TRUE, 
  sims = 1000)

# Step 4) Summarize mediation results (Appendix E, Table E11)
summary(medmodel.trait.nnr) # attribute: personality/leadership traits (Model C7.1)
summary(medmodel.settg.nnr) # attribute: incivility setting (Model C7.2)

